Code Example

Prepare Environment

Load necessary modules and attempt to compile the GrOpt library if it isn't already

%load_ext watermark

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import chart_studio.plotly as py
import cufflinks as cf
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go

# This attempts to re-compile the library in case it has been changed, mostly for debug, but won't do anything
# if nothing is changed
import build_gropt
build_gropt.build_gropt()
import gropt

from helper_utils import *
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from IPython.core.display import display, HTML
init_notebook_mode(connected = True)
config={'showLink': False, 'displayModeBar': False}
from timeit import default_timer as timer

%matplotlib inline

Building GrOpt . . .

Print all used libraries and system info, so as to enable reproducibility

%watermark -v -m -p watermark,pandas,numpy,matplotlib,chart_studio,cufflinks,seaborn,plotly
print()
%watermark -u -n -t -z

Python implementation: CPython
Python version       : 3.7.8
IPython version      : 7.16.1

watermark   : 2.1.0
pandas      : 1.0.1
numpy       : 1.18.1
matplotlib  : 3.1.3
chart_studio: 1.1.0
cufflinks   : 0.17.3
seaborn     : 0.10.0
plotly      : 4.14.3

Compiler    : GCC 7.5.0
OS          : Linux
Release     : 4.19.150+
Machine     : x86_64
Processor   : x86_64
CPU cores   : 8
Architecture: 64bit


Last updated: Sat Jan 30 2021 00:23:39UTC

Interactive figures

First define the function that will plot waveforms interactively

def plot_waveform_interactive(G, params, plot_moments = True, plot_eddy = True, plot_pns = True, plot_slew = True,
                  suptitle = '', eddy_lines=[], eddy_range = [1e-3,120,1000]):
    
    n_axis = params.get('Naxis', 1)
    num_traces = n_axis + int(plot_slew)*n_axis + int(plot_moments)*3 + int(plot_pns) + int(plot_eddy)
    current_trace = 0

    TE = params['TE']
    T_readout = params['T_readout']
    diffmode = 0
    if params['mode'][:4] == 'diff':
        diffmode = 1
    
    dt = (TE-T_readout) * 1.0e-3 / G.shape[1]
    tt = np.arange(G.shape[1]) * dt * 1e3
    tINV = TE/2.0
    
    bval = get_bval(G, params)
    blabel = 'b-value=%.0f' % bval # TODO: should add "$mm^{2}/s$"
    
    fig_title = None
    if suptitle:
        fig_title = suptitle + ': ' + blabel
    elif diffmode > 0:
        fig_title = blabel
        
    fig = go.Figure(layout=dict(xaxis=dict(title="$t [ms]$")))
    buttons = []
    
    # Gradient
    for i in range(n_axis):
        fig.add_trace(go.Scatter(x=tt, y=G[i]*1000, showlegend=False))
    buttons.append(dict(label = "Gradient",
                   method = "update",
                   args = [{"visible": [True if x>= current_trace and x < current_trace + n_axis
                                        else False for x in range(num_traces)],
                           "showlegend":False},
                           {"shapes":[],
                           "xaxis": {"title": "$t [ms]$"}}]))
    current_trace += n_axis
    
    # Slew
    if plot_slew:
        for i in range(n_axis):
            fig.add_trace(go.Scatter(x=tt[:-1], y=np.diff(G[i])/dt, visible=False, showlegend=False))
        buttons.append(dict(label = "Slew",
                       method = "update",
                       args = [{"visible": [True if x>= current_trace and x < current_trace + n_axis
                                            else False for x in range(num_traces)],
                               "showlegend":False},
                               {"shapes":[],
                               "xaxis": {"title": "$t [ms]$"}}]))
        current_trace += n_axis
            
    # Moments
    if plot_moments:
        moment_lines = []
        moment_lines.append({"type":"line", "x0":tt[0], "x1":tt[-1], "y0":0, "y1":0,
                          "line":{"color":"#777777", "width":1, "dash":"dash"}})
        mm = get_moment_plots(G, T_readout, dt, diffmode)
        for i in range(3):
            if diffmode == 1:
                mmt = mm[i]/np.abs(mm[i]).max()
            if diffmode == 0:   
                if i == 0:
                    mmt = mm[i]*1e6
                if i == 1:
                    mmt = mm[i]*1e9
                if i == 2:
                    mmt = mm[i]*1e12
            fig.add_trace(go.Scatter(x=tt, y=mmt, visible=False, name='$M_{%d}$'%i))
        buttons.append(dict(label = "Moments",
                       method = "update",
                       args = [{"visible": [True if x>= current_trace and x < current_trace + 3
                                            else False for x in range(num_traces)],
                               "showlegend":True},
                               {"shapes":moment_lines,
                               "xaxis": {"title": "$Time [ms]$"}}]))
        current_trace += 3
            
    # Eddy
    if plot_eddy:
        all_lam = np.linspace(eddy_range[0],eddy_range[1],eddy_range[2])
        all_e = []
        for lam in all_lam:
            lam = lam * 1.0e-3
            r = np.diff(np.exp(-np.arange(G[0].size+1)*dt/lam))[::-1] # TODO: 3-axis case, right now just assumes 1 axis
            all_e.append(100*r@G[0])
        fig.add_trace(go.Scatter(x=all_lam, y=all_e, visible=False, showlegend=False))

        eddy_draw = []
        for e in eddy_lines:
            min_e = min(all_e)
            max_e = max(all_e)
            amp = 0.1*max(abs(min_e), abs(max_e))
            eddy_draw.append(dict(type="line",
                x0=e, x1=e, y0=min(all_e)-amp, y1=max(all_e)+amp,
                line=dict(color="#D64B4B",
                    width=1,
                    dash="dot")))
        eddy_draw.append({"type":"line", "x0":all_lam[0], "x1":all_lam[-1], "y0":0, "y1":0,
                          "line":{"color":"#777777", "width":1, "dash":"dash"}})
        buttons.append(dict(label = "Eddy",
                       method = "update",
                       args = [{"visible": [True if x>= current_trace and x < current_trace + 1
                                            else False for x in range(num_traces)],
                               "showlegend":False},
                               {"shapes":eddy_draw,
                               'xaxis': {'title': "$\lambda [ms]$"}}]))
        current_trace += 1
    
    # PNS
    if plot_pns:
        pns_lines = []
        pns_lines.append({"type":"line", "x0":tt[0], "x1":tt[-2], "y0":1, "y1":1,
                      "line":{"color":"#D64B4B", "width":1, "dash":"dot"}})
        pns_lines.append({"type":"line", "x0":tt[0], "x1":tt[-2], "y0":0, "y1":0,
                          "line":{"color":"#777777", "width":1, "dash":"dash"}})
        pns = np.abs(get_stim(G, dt))
        fig.add_trace(go.Scatter(x=tt[:-1], y=pns, visible=False, showlegend=False))
        buttons.append(dict(label = "PNS",
                       method = "update",
                       args = [{"visible": [True if x>= current_trace and x < current_trace + 1
                                            else False for x in range(num_traces)],
                               "showlegend":False},
                               {"shapes":pns_lines,
                               "xaxis": {"title": "$Time [ms]$"}}]))
        current_trace += 1
        
    # Layout
    update_menu = [
        dict(active = 0,
            direction="down",
            pad={"r": 10, "t": 10},
            showactive=True,
            x=1.025,
            xanchor="left",
            y=.85,
            yanchor="top",
            buttons=buttons
            )]
    
    fig.update_layout(title=fig_title,
        updatemenus=update_menu,
        width=900,
        height=550,
        autosize=False,
        margin=dict(t=50, b=50, l=0, r=50),
        template="plotly_white")
    
    plot(fig, filename = 'fig.html', config = config)
    display(HTML('fig.html'))

Test the function

params = {}
params['mode'] = 'free'
params['gmax']  = 0.05
params['smax']  = 50.0
params['moment_params']  = [[0, 0, 0, -1, -1, 11.74, 1.0e-3]]
params['moment_params'].append([1, 0, 0, -1, -1, -11.74, 1.0e-3])
params['moment_params'].append([2, 0, 0, -1, -1, 11.74, 1.0e-3])
params['TE']  = 1.0
params['dt']  = 20e-6
params['Naxis'] = 3

G, dd = gropt.gropt(params, verbose=1)

plot_waveform_interactive(G, params)

params = {}
# Maximize b-value for diffusion waveforms
params['mode'] = 'diff_bval'

# Hardware constraints
params['gmax']  = 50.0 # Max Gradient Amplitude [mT/m], you can use T/m here too
params['smax']  = 50.0 # Max Slewrate [mT/m/ms]

# Moment nulling
params['MMT']  = 0

# Sequence TE and dt of output [ms]
params['TE']  = 60
params['dt']  = 400e-6

# Time from end of diffusion waveform to TE [ms]
params['T_readout']  = 16.0
# Time of excitation 90 [ms]
params['T_90']  = 4.0
# Time for 180 flip [ms]
params['T_180']  = 6.0

# Run optimization
G, dd = gropt.gropt(params, verbose=1)

plot_waveform_interactive(G, params, eddy_lines=[50])

params = {}
params['mode'] = 'free'
params['gmax']  = 0.08
params['smax']  = 200.0
params['moment_params']  = [[0, 0, 0, -1, -1, 11.74, 1.0e-3]]
params['moment_params'].append([1, 0, 0, -1, -1, -11.74, 1.0e-3])
params['moment_params'].append([2, 0, 0, -1, -1, 11.74, 1.0e-3])
params['TE']  = 0.675
params['dt']  = 10e-6
params['Naxis'] = 3
params['pns_thresh'] = 1.0

G, dd = gropt.gropt(params, verbose=1)

plot_waveform_interactive(G, params, plot_slew=False)

params = {}
params['mode'] = 'diff_bval'
params['gmax']  = 0.08
params['smax']  = 100.0
params['MMT']  = 0
params['TE']  = 50.0
params['T_readout']  = 16.0
params['T_90']  = 3
params['T_180']  = 5
params['dt']  = 200e-6

G, dd = gropt.gropt(params, verbose=1)

plot_waveform_interactive(G, params)

params = {}
params['mode'] = 'diff_bval'
params['gmax']  = 0.05
params['smax']  = 50.0
params['MMT']  = 2
params['TE']  = 97.0
params['T_readout']  = 16.0
params['T_90']  = 4.0
params['T_180']  = 6.0
params['dt']  = 200e-6

G, dd = gropt.gropt(params, verbose=1)

plot_waveform_interactive(G, params, eddy_lines=[50])

params = {}
params['mode'] = 'free' # Free mode indicates we are in a feasibility search, i.e. no objective function
params['dt'] = 4e-6     # Raster time of the gradient waveform being optimized
params['gmax'] = 50     # Maximum gradient amplitude, mT/m
params['smax'] = 120.0  # Maximum slew rate, mT/m/ms 
params['moment_params'] = [[0, 0, 0, -1, -1, 11.75, 1.0e-3]] # Constraining for M0 = 11.75 with a tolerance of 1.0e-3


G, T_min = get_min_TE(params, max_TE = 2, verbose = 1)
print('Waveform duration =', round(T_min,2), 'ms')
plot_waveform_interactive(G, params, plot_moments = True, plot_eddy = False, plot_pns = False, plot_slew = True)

Testing TE = 1.050 0.575 0.812 0.694 0.634 0.605 0.620 0.627 0.631 0.629 0.630 Final TE = 0.631 ms
Waveform duration = 0.63 ms